gisR

GIS-Review

Thomas Knecht, Ueli Mauch

11. Juli 2024

Was ist R

R ist eine offene Programmiersprache welche ursprünglich für statistische Zwecke erstellt wurde.

Der Einsatzbereich von R hat sich heutzutage weit über die statistischen Anwendungen etabliert.

Hinweis

Gerade für Data-Wrangling und Data-Engineering kommt die Stärke von R besonders zu tragen.

GIS mit R

Seit einiger Zeit hat sich R im Bereich GIS weiterentwickelt. Es lässt sich nun sehr einfach GIS-Prozesse mit Datenanalyse verbinden und es ist kein Medienbruch mehr nötig.

Packages (gleich wie libraries in python)

  • sf: Arbeiten mit Feature Classes

  • terra: Arbeiten mit Raster-Daten

  • stars: Arbeiten mit Spatio-Temporal Arrays

Einlesen der Daten

# Lesen der Gebäudedaten
gwr <- read.csv(paste0("https://www.web.statistik.zh.ch/ogd/daten/",
                       "ressourcen/KTZH_00002022_00004064.csv")
                ) |>
  dplyr::rename_all(tolower)


# Umwandeln in ein sf-Objekt -> Geodatenobjekt
gwr_sf <- sf::st_as_sf(
  gwr, 
  coords = c("e.gebaeudekoordinate", "n.gebaeudekoordinate"), 
  crs = 2056
)

# interaktives Darstellen
mapview::mapview(head(gwr_sf, 100))

gemeinden <- sf::st_read(
  paste0("https://maps.zh.ch/wfs/OGDZHWFS?SERVICE=WFS&VERSION=2.0.0&",
         "Request=getfeature&TYPENAME=ms:ogd-0095_arv_basis_up_gemeinden_seen_f&",
         "outputformat=geojson")
  ) |>
  dplyr::rename_all(tolower)

Analysieren der Daten

# Matchen der Gemeindeinformation an die GWR-Daten
# Die Geometrie wird hier entfernt, da die Datenanalyse danach 
# viel schneller läuft
gwr_sf_gem <- gwr_sf |>
  dplyr::select(
    gebaeudekategorie_code, 
    gebaeudekategorie_bezeichnung
    ) |>
  sf::st_join(gemeinden) |>
  sf::st_drop_geometry()


# Berechnung der Anteile der Gebäude ohne Wohnnutzung 
# am gesamten Gebäudebestand einer Gemeinde
geb_ohne_wohnnutzung_pro_gem <- gwr_sf_gem |>
  dplyr::group_by(
    gebaeudekategorie_code, 
    gebaeudekategorie_bezeichnung, 
    bfs, 
    gemeindename
  ) |>
  dplyr::summarise(anzahl = dplyr::n()) |>
  dplyr::ungroup() |>
  dplyr::group_by(bfs, gemeindename) |>
  dplyr::mutate(anteil = round(anzahl/sum(anzahl)*100, 2)) |>
  dplyr::ungroup() |>
  dplyr::filter(gebaeudekategorie_code == 1060)

Vorbereiten für die Karten

# Hinzufügen der Gemeindepolygone
# Herausfiltern der Seen und fixen der Polygone
geb_ohne_wohnutzung_sf <- gemeinden  |>
  dplyr::select(bfs) |>
  dplyr::left_join(geb_ohne_wohnnutzung_pro_gem, by = "bfs") |>
  dplyr::filter(bfs != 0) |>
  sf::st_make_valid() |>
  # tmap vewendet die erste Spalte als Bezeichner, 
  #deshalb hier die Änderung in der Spaltenanordnung
  dplyr::select(gemeindename, dplyr::everything())


geb_ohne_wohnnutzung_sf_generalized <- geb_ohne_wohnutzung_sf |>
  rmapshaper::ms_simplify(keep = 0.005)

Statische Karte

## Nur für Präsentation!
geb_ohne_wohnnutzung_sf_generalized <- readRDS(here::here("extdata/data.RDS"))

# Statisches Darstellung
tmap::tm_shape(geb_ohne_wohnnutzung_sf_generalized) +
  tmap::tm_polygons("anteil") +
  tmap::tm_compass() +
  tmap::tm_scale_bar()

Interaktive Karte

## Nur für Präsentation!
geb_ohne_wohnnutzung_sf_generalized <- readRDS(here::here("extdata/data.RDS"))

# Setzt tmap auf interaktiv
tmap::tmap_mode("view")

# Interaktive Darstellung
tmap::tm_shape(geb_ohne_wohnnutzung_sf_generalized) +
    tmap::tm_polygons("anteil")

Remote Sensing mit R

Extensiv genutzte Wiesen müssen mindestens einmal pro Jahr gemäht werden und das Schnittgut muss abgeführt werden. Die Flächen dürfen in Abhängigkeit der Zone jeweils frühestens Mitte Juni bis Mitte Juli genutzt werden1

Fragestellung

Gibt es Landwirtschaftsflächen1 mit Schnittzeitpunkt 15. Juni, welche bereits vorher gemäht wurden?

Workflow:2

  • Daten herunterladen
  • NDVI berechnen
  • Statistik für Landwirtschaftsflächen errechnen

STAC

  • SpatioTemporal Asset Catalog1

  • Wir benutzen den Planetary Computer Data Catalog von Microsoft2

STAC

  • Collections abfragen
stac_source <- rstac::stac(
  "https://planetarycomputer.microsoft.com/api/stac/v1",
  force_version = "1.0.0"
)
collections_query <- stac_source |>
  rstac::collections()

rstac::get_request(collections_query)
###Collections
- collections (123 item(s)):
  - daymet-annual-pr
  - daymet-daily-hi
  - 3dep-seamless
  - 3dep-lidar-dsm
  - fia
  - sentinel-1-rtc
  - gridmet
  - daymet-annual-na
  - daymet-monthly-na
  - daymet-annual-hi
  - ... with 113 more collection(s).
- field(s): collections, links

STAC

# BBOX definieren
ktzh_bbox <- sf::read_sf("geodata/Gemeindegrenzen_-OGD.gpkg") |> 
  sf::st_transform("EPSG:4326") |> 
  sf::st_make_valid() |> 
  sf::st_union() |> 
  sf::st_bbox() |> 
  format_bbox()
  

# STAC Query bauen
stac_query <- rstac::stac_search(
  q = stac_source,
  collections = "sentinel-2-l2a",
  datetime = "2023-06-10/2023-06-14",
  bbox = ktzh_bbox
) |>
  rstac::ext_filter(
    `eo:cloud_cover` <= 20)

Filtern nach Gebiet, Zeitraum, Collection Sentinel-2-L2A. Bewölkung einschränken.

Mehr zum Download in scripts/dlFromSTAC.R

Geodaten von WFS laden

wfs_zh <- "https://maps.zh.ch/wfs/OGDZHWFS"

url <- httr::parse_url(wfs_zh)
url$query <- list(service = "WFS",
                  version = "2.0.0",
                  request = "GetFeature",
                  typename = "ms:ogd-0095_arv_basis_up_bezirke_f",
                  outputformat = "geojson")

request <- httr::build_url(url)
bezirke <- sf::read_sf(request)

# Save Bezirke Layer
sf::write_sf(bezirke, "geodata/Bezirke.gpkg")

WFS GetFeature Request, am Beispiel der OGD Bezirksgrenzen

Raster

  • Bänder-Kombination & Visualisierung (terra Package)
RGB FCIR
RGB FCIR

NDVI

Quelle

Rasterrechner

  • 4-Band Raster als Basis
  • NDVI nach dieser Formel

\[ NDVI = \frac{NIR - R}{NIR + R} \]

  • Berechnung simpel:
# calculate NDVI using the red (band 1) and nir (band 4) bands
sent2_ndvi <- (RGBI_zh[[4]] - RGBI_zh[[1]]) / (RGBI_zh[[4]] + RGBI_zh[[1]])

sent2_ndvi |> terra::writeRaster("geodata/ndvi_terra.tif")

Statistik

  • Nutzflächen auswählen, Feld area berechnen
  • 2 Filter anwenden
# Load Polygon Layer
LWNutz_Diels <- sf::read_sf("geodata/LWNutz_Dielsdorf.gpkg") |>
  dplyr::mutate(area = sf::st_area(geom))

# Set aller Flächen, die erst am 15.06. sollten geschnitten werden:
BF_Fl <- LWNutz_Diels |>
  dplyr::filter(harvest_date == "15.06.") |>
  # Drop die kleinsten Flächen unter 500 m^2
  dplyr::filter(area > units::set_units(500, "m^2"))

Zonale Statistik

  • Statistik pro Nutzfläche, jeweils der mean wird berechnet
bfs_parz area mean_ndvi geom
0088_1041 2917.503 [m^2] 0.5237937 POLYGON ((2678319 1262664, ...
0088_965 10950.015 [m^2] 0.5349045 POLYGON ((2678276 1262340, ...
0100_1517 4092.843 [m^2] 0.5283303 POLYGON ((2678039 1265958, ...
0088_951 1817.895 [m^2] 0.5004897 POLYGON ((2678444 1262461, ...

Verteilung

Map

# Setzt tmap auf interaktiv
border <- c(2669379, 1252226, 2684597, 1269716)
breaks <- c(-0.4, -0.2, 0, 0.2, 0.4, 1)

tmap::tm_shape(BF_Fl, bbox = border) +
tmap::tm_polygons("mean_ndvi", style = "fixed", breaks = breaks) +
tmap::tm_shape(BF_Fl |> dplyr::filter(mean_ndvi <= 0.4), bbox = border, name = "NDVI <= 0.4 ") +
tmap::tm_polygons("mean_ndvi", style = "fixed", breaks = breaks, legend.show = FALSE)
BF_Fl |> dplyr::filter(`mean_ndvi` <= 0.4) |> 
  dplyr::select(blw_name, region, mean_ndvi, area) |>
  sf::st_drop_geometry() |> head(10) |> 
  knitr::kable() |> kableExtra::kable_styling(font_size = 22)
blw_name region mean_ndvi area
Extensiv genutzte Wiesen (ohne Weiden) Chrüzstrasse 0.3896368 578.9319 [m^2]
Hecken, Feld-, Ufergehölze mit Krautsaum Riedthof Ost 0.3903728 798.4844 [m^2]
Wenig intensiv genutzte Wiesen (ohne Weiden) Hard Humusdepot 0.3058146 10962.5786 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Lierehölzli 0.3605108 681.2977 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Fallmet 0.3872881 2092.0953 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Hodleten 0.3928246 1726.3641 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Spitz 0.3838461 1914.2668 [m^2]
Hecken, Feld-, Ufergehölze mit Krautsaum Heitlig Steinhaufen 0.3893114 917.3088 [m^2]
Hecken, Feld-, Ufergehölze mit Krautsaum Talacher 0.3867528 1794.5135 [m^2]
Extensiv genutzte Wiesen (ohne Weiden) Geeren 0.3962817 872.1179 [m^2]

Diskussion